individual_level = 
  individual_level %>%
  mutate(sc02_14_14 = str_pad(sc02_14_14, width = 2, side = "left", pad = "0")) %>%
  mutate(kab_code = as.numeric(paste0(sc01_14_14, sc02_14_14)))

crosswalk = read.csv("./_3_data/crosswalks/district_crosswalk.csv") %>% dplyr::select(bps_1993, bps_2014)

individual_level <- left_join(individual_level, crosswalk, by = c("kab_code" = "bps_2014"))

#2014 District Fixed Effects

band <- individual_level %>% filter(abs(forcing) < 10000)
model_df1_1 <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing + factor(kab_code), data=band)
model_df1_2 <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing + factor(kab_code), data=band)
model_df1_3 <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing + factor(kab_code), data=band)
model_df1_4 <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing + factor(kab_code), data=band)

model_df1_1 <- coeftest(model_df1_1, vcov=cluster.vcov(model_df1_1, cluster = band$mkid14))
model_df1_2 <- coeftest(model_df1_2, vcov=cluster.vcov(model_df1_2, cluster = band$mkid14))
model_df1_3 <- coeftest(model_df1_3, vcov=cluster.vcov(model_df1_3, cluster = band$mkid14))
model_df1_4 <- coeftest(model_df1_4, vcov=cluster.vcov(model_df1_4, cluster = band$mkid14))

model_df1_1 <- tidy(model_df1_1) %>% mutate(outcome = "outcome1", level = "2014 District FE")
model_df1_2 <- tidy(model_df1_2) %>% mutate(outcome = "outcome2", level = "2014 District FE")
model_df1_3 <- tidy(model_df1_3) %>% mutate(outcome = "outcome3", level = "2014 District FE")
model_df1_4 <- tidy(model_df1_4) %>% mutate(outcome = "outcome4", level = "2014 District FE")


#1993 District Fixed Effects

band <- individual_level %>% filter(abs(forcing) < 10000)

model_df2_1 <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing + factor(bps_1993), data=band)
model_df2_2 <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing + factor(bps_1993), data=band)
model_df2_3 <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing + factor(bps_1993), data=band)
model_df2_4 <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing + factor(bps_1993), data=band)

model_df2_1 <- coeftest(model_df2_1, vcov=cluster.vcov(model_df2_1, cluster = band$mkid14))
model_df2_2 <- coeftest(model_df2_2, vcov=cluster.vcov(model_df2_2, cluster = band$mkid14))
model_df2_3 <- coeftest(model_df2_3, vcov=cluster.vcov(model_df2_3, cluster = band$mkid14))
model_df2_4 <- coeftest(model_df2_4, vcov=cluster.vcov(model_df2_4, cluster = band$mkid14))

model_df2_1 <- tidy(model_df2_1) %>% mutate(outcome = "outcome1", level = "1993 District FE")
model_df2_2 <- tidy(model_df2_2) %>% mutate(outcome = "outcome2", level = "1993 District FE")
model_df2_3 <- tidy(model_df2_3) %>% mutate(outcome = "outcome3", level = "1993 District FE")
model_df2_4 <- tidy(model_df2_4) %>% mutate(outcome = "outcome4", level = "1993 District FE")

#2014 provincial fixed effects
band <- individual_level %>% filter(abs(forcing) < 10000)
model_df3_1 <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing + factor(sc01_14_14), data=band)
model_df3_2 <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing + factor(sc01_14_14), data=band)
model_df3_3 <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing + factor(sc01_14_14), data=band)
model_df3_4 <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing + factor(sc01_14_14), data=band)

model_df3_1 <- coeftest(model_df3_1, vcov=cluster.vcov(model_df3_1, cluster = band$mkid14))
model_df3_2 <- coeftest(model_df3_2, vcov=cluster.vcov(model_df3_2, cluster = band$mkid14))
model_df3_3 <- coeftest(model_df3_3, vcov=cluster.vcov(model_df3_3, cluster = band$mkid14))
model_df3_4 <- coeftest(model_df3_4, vcov=cluster.vcov(model_df3_4, cluster = band$mkid14))

model_df3_1 <- tidy(model_df3_1) %>% mutate(outcome = "outcome1", level = "2014 Province FE")
model_df3_2 <- tidy(model_df3_2) %>% mutate(outcome = "outcome2", level = "2014 Province FE")
model_df3_3 <- tidy(model_df3_3) %>% mutate(outcome = "outcome3", level = "2014 Province FE")
model_df3_4 <- tidy(model_df3_4) %>% mutate(outcome = "outcome4", level = "2014 Province FE")


#make plot
plot <-
  bind_rows(model_df1_1, model_df1_2, model_df1_3, model_df1_4,
          model_df2_1, model_df2_2, model_df2_3, model_df2_4,
          model_df3_1, model_df3_2, model_df3_3, model_df3_4) %>%
  mutate(outcome_label = case_when(outcome == "outcome1" ~ "Place of worship",
                                   outcome == "outcome2" ~ "Live in village",
                                   outcome == "outcome3" ~ "Live in neighborhood",
                                   outcome == "outcome4" ~ "Rent room",
                                   TRUE ~ NA_character_)) %>%
  filter(term == "opium_legal") %>%
  #mutate(Inequality = level) %>%
  ggplot(aes(x=estimate, y = outcome_label, group = level, shape = level, color = level)) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbarh(aes(xmin=estimate-1.67*std.error, xmax = estimate+1.67*std.error), position = position_dodge(width = 0.4), height = 0) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  xlab("Effect of Opium Concession System") +
  theme_bw() +
  scale_color_grey() +
  theme(axis.ticks.y.left = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title.y = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.title = element_blank(),
        axis.line.y.left = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))


ggsave("./_4_outputs/figures/appendix_4_geo_fe_plot.png", plot = plot,  width = 7, height = 4)

